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ABSTRACT 

Magnetic interactions between a protostar and its accretion disc can induce warping 
in the disc and produce secular changes in the stellar spin direction, so that the spin 
axis may not always be perpendicular to the disc. This may help explain the 7-degree 
misalignment between the ecliptic plane of the solar system and the sun's equatorial 
planem as well as play a role in producing the recently observed spin-orbit misalign- 
ment in a number of exoplanetary systems. We study the dynamics of warped pro- 
toplanetary discs under the combined effects of magnetic warping/precession torques 
and internal stresses in the disc, including viscous damping of warps and propaga- 
tion of bending waves. We show that when the outer disc axis is misaligned with the 
stellar spin axis, the disc evolves towards a warped steady-state on a timescale that 
depends on the disc viscosity or the bending wave propagation speed, but in all cases 
is much shorter than the timescale for the spin evolution (of order of a million years) . 
Moreover, for the most likely physical parameters characterizing magnetic protostars, 
circumstellar discs and their interactions, the steady-state disc, averaged over the stel- 
lar rotation period, has a rather small warp such that the whole disc lies approximately 
in a single plane determined by the outer disc boundary conditions, although more 
extreme parameters may give rise to larger disc warps. In agreement with our recent 
analysis ( |Lai et aL||2010 ) based on flat discs, we find that the back-reaction magnetic 



torques of the slightly warped disc on the star can either align the stellar spin axis with 
the disc axis or push it towards misalignment, depending on the parameters of the 
star-disc system. This implies that newly formed planetary systems may have a range 
of inclination angles between the stellar spin axis and the orbital angular momentum 
axis of the planetary orbits. 

Key words: accretion, accretion discs - planetary systems: protoplanetary discs - 
stars: magnetic fields 



1 INTRODUCTION 



In a recent paper Lai et al. (2010), hereafter Paper I], we 



proposed a novel mechanism for producing misalignment be- 
tween the spin axis of a protostar and the normal vector of 
its circumstellar disc. Our work was motivated by recent 
measurements of the sky-projected stellar obliquity using 
the Rossiter-McLaughlin effect in transiting exoplanetary 
systems, which showed that a large fraction of the systems 
containing hot Jupiters have misaligned stellar spin with re- 
spect to the planetary angular momentum axis [see |Triaud| 
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et al. (20101; Winn et al. (20101 and references therein] 



Additional evidence for nonzero stellar obliquity came from 
the statistical analysis of the apparent rotational velocities 
(usinu) of planet-bearing stars ( |Schlaufmanl[2uTu) . 



The basic mechanism ("Magnetically driven misalign- 
ment") for producing spin - disc misalignment in accreting 
protostellar systems can be sumarized as follows (Paper I). 
The magnetic field of a protostar (with B+ <; 10 3 G) pene- 
trates the inner region of its accretion disc. These field lines 
link the star and the disc in a quasi-cyclic fashion (e.g., mag- 
netic field inflation followed by reconnection; see |Bouvier et| 
[af] ( |2007| > and |Alencar et al.| ( |2010[ | for observational evi- 
dence). Differential rotation between the star and the disc 
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not only leads to the usual magnetic braking torque on the 
disc, but also a warping torque which tends to push the 
normal axis of the inner disc away from the spin axi^] Hy- 
drodynamical stresses in the disc, on the other hand, tend 
to inhibit significant disc warping. The result is that, for a 
given disc orientation imposed at large radii (e.g., by the an- 
gular momentum of the accreting gas falling onto the disc), 
the back-reaction of the warping torque can push the stel- 
lar spin axis toward misalignment with respect to the disc 
normal vector. Planets formed in the disc will then have a 
misaligned orbital normal axis relative to the stellar spin 
axis, assuming that no evolution mechanism occurring af- 
ter the dissipation of the disc forces the alignment of the 
system. 

The process of planetary system formation can be 



roughly divided into two stages (Juric & Tremaine 20081 



In the first stage, which lasts a few million years until the 
dissipation of the gaseous protoplanetary disc, planets are 
formed and undergo migration due to tidal interactions with 
the gaseous disc Lin et al. (19961; see Papaloizou et al. 



(20071 for a review]. The second stage, which lasts from 



when the disc has dissipated to the present, involves dy- 
namical gravitational interactions between multiple planets, 
if they are produced in the first stage in a sufficiently close- 
packed configuration (Juric & Tremaine 2008; Chatterjee et 



aL||2008 1, and/or secular interactions with a distant planet 



or stellar companion (Eggleton & Kiseleva-Eggleton 



2001 



Wu fc Murray|2003l |Fabrycky fc Tremaine||2007| |Wu et al.| 



2007| Nagasawa et al.|2008 ) . The eccentricity distribution of 
exoplanetary systems and the recent observational results 
on the spin - orbit misalignment suggest that the physi- 
cal processes in the second stage play an important role in 
determining the properties of exoplanetary systems. Never- 
theless, the importance of the first stage cannot be neglected 
as it sets the initial condition for the possible evolution in 
the second stage. Our result in Paper I shows that at the 
end of the first stage, the symmetry axis of the planetary 
orbit may be inclined with respect to the stellar spin axis. 

At first sight, it may seem strange that the magnetic 
field effects can drive the stellar spin axis toward misalign- 
ment with respect to the disc symmetry axis, given that the 
spin angular momentum of the star ultimately comes from 
the disc and the disc contains a large reservoir of angular 
momentum. The key to understand this is to realize that 
when the gas reaches the magnetosphere boundary, its an- 
gular momentum is much smaller than in the outer disc (the 
specific angular momentum of the disc is jdisc(f*) = V GMr 
for a Keplerian disc), and any magnetic torque, which in 
general can break the axisymmetry of the system, is of the 
same order of magnitude as the accretion torque on the star. 

A key assumption adopted in Paper I for the calculation 
of the magnetic torque on the star from the disc is that the 
disc is flat. This is a nontrivial assumption. Indeed, the mag- 
netic coupling between the star and the disc operates only 
in the innermost disc region (e.g., between the inner radius 
rin and n n t ~ 1.5n n ), and this region has a much smaller 



1 The warping torque vanishes if the angular momentum of the 
disc is exactly aligned with the stellar spin, but exists for arbi- 
trary small angles. A flat disc in the aligned configuration is in 
an unstable equilibrium 



moment of inertia than the star. Therefore, if there were no 
coupling between this inner disc region and the outer disc, 
the inner disc would be significantly warped on a timescale 
much shorter than the timescale for changing the stellar spin 



(Pfeiffer & Lai 2004h. If there is any secular change in the 



stellar spin direction, the inner disc warp would then fol- 
low the varying spin axis. Clearly, in order to determine the 
long-term spin evolution of the star, it is important to under- 
stand the dynamics of the warped disc, taking into account 
the magnetic torques on the inner disc and the hydrody- 
namical coupling between different disc regions. This is the 
goal of our paper. 

To be more specific, there is a hierarchy of timescales 
related to the combined evolution of the stellar spin and the 
disc warp: 

(i) The dynamical time tdyn associated with the spin fre- 
quency U) s , disc rotation frequency f2 and the beat frequency 
cj s — Q\. This is much shorter than the effects (steady-state 
disc warping and spin evolution) we study in this paper. 

(ii) The warping/precession timescale of the inner disc 
[see Eq. @] 
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where M*, R*, -B* are the mass, radius and surface (dipole) 
magnetic field of the protostar, respectively, (9* is the incli- 
nation angle of the stellar dipole relative to the spin, E is 
the disc surface density, and C, is a dimensionless magnetic 
twist parameter of order unity related to the strength of the 
azimuthal magnetic field generated by star-disc twist. 

(iii) The disc warp evolution timescale £disc • This is the 
time for the disc to reach a steady-state under the com- 
bined effects of magnetic torques and internal fluid stresses 
(see Section 5). For high-viscosity discs, tdisc is the viscous 
diffusion time for the disc warp [see Eq. ( 34 1] and depends 
on the viscosity parameter a and the disc thickness S = H/r: 



; vis ~(3000yrs)(^)(A) "(_!_) 



3/2 



(2) 



For low- viscosity discs (a & 5), idisc is the propagation time 
of bending waves across the whole disc and depends on the 
sound speed. In general, idisc can be several orders of mag- 
nitude larger than t w . 

(iv) The stellar spin evolution timescale. The magnetic 
misalignment torque on the star is of order p 2 /rf n (fj, being 
the magnetic dipole moment of the star) , which is compara- 
ble to the fiducial accretion torque, given for Keplerian discs 
by J\fo = M^GM+r^. Assuming the spin angular momen- 
tum J s = Q.2M*R*U] S (the value for a V = 5/3 polytrope, 
representing a convective star), we find the spin evolution 
time 



^(1.25Myr)(^- 



M 



lO- 8 M yr- 



(3) 



In general i sp in 2> idisc- In this paper we will study the 
evolution of the disc warp on timescales ranging from t w 
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to tdisc, and the evolution of the stellar spin direction on 
timescales of order i ap in- 

It is important to note that we are not interested in disc 
warpings that vary on the dynamical timescale tdyn in this 
paper. In general, when the stellar dipole axis is inclined 
with respect to the spin axis, there will be periodic verti- 
cal forces at the rotation frequency of the star acting on 
the inner discQ These periodic forces will lead to the warp- 
ing of the disc, particularly for low-viscosity discs in which 
bending waves propagate ( |Terquem fc Papaloizou|2000| |Lai| 
& Zhang 2008). Indeed, there is observational evidence for 
such magnetically-warped discs. For example, the recurrent 
luminosity dips observed in the classical T Tauri star AA 
Tauri has been attributed to the periodic occultation of the 



central star by a warped inner disc (Bouvier et al. 20071 



However, such dynamical disc warps average exactly to zero 
over a rotation period and have no effect on the secular evo- 
lution of the system. 

The remainder of the paper is organized as follows. In 
Section |2j we summarize our analytical model of magne- 
topshere - disc interaction and derive the equation for the 
evolution of the stellar spin axis when the disc is warped. 
In Section 3 we present theoretical formalisms for deter- 
mining the steady-state and time evolution of warped discs, 
for both high-viscosity regime (where warps propagate dif- 
fusively) and low-viscosity regime (where warps propagate 
as bending waves). An approximate analytical expression 
for the steady-state linear warp is also derived (see Sec- 
tion 3.2.2). In Section 4 we present numerical results for 
the steady-state disc warp profiles under various conditions 
and in Section 5 we study the time evolution of disc warps. 
We examine in Section 6 how the inner disc warp and the 
stellar spin evolution respond to variations of the outer disc, 
and discuss in Section 7 how this could, in principle, lead to 
anti-aligned planetary orbits, even for discs with initial an- 
gular momentum nearly aligned with the stellar spin. We 
conclude in Section 8 with a discussion of our results. 



2 ANALYTIC MODEL OF THE DISC - 
MAGNETIC STAR SYSTEM 

2.1 Magnetic Torques on the Disc 

The interaction between a magnetic star and a disc is com- 
plex (see references in Paper I). However, the key physical 
effects of this interaction on the disc can be described ro- 
bustly in a parametrized manner. The model used through- 
out this paper is detailed in Paper I. Here, we will limit 
ourselves to a brief summary of the magnetic torques acting 
on the disc. 

The stellar magnetic field disrupts the accretion disc 
at the magnetospheric boundary, where the magnetic and 
plasma stresses balance. For a dipolar magnetic field with 
magnetic moment fi, we have 



fi. 



1/7 



(4) 



2 The forcing frequency may also be twice of the spin frequency 
under certain conditions (e.g., when the disc is partially diamag- 



where r\ is a dimensionless constant somewhat less than 
unity (77 ~ 0.5 according to recent numerical simulations; 
see Long et al. 2005^]). We take r- m to be the inner edge of 
the disc. Before being disrupted, the disc generally experi- 
ences nontrivial magnetic torques from the star (Lai 1999; 
Paper I). Consider a cylindrical coordinate system (r,(f>,z), 
with the vertical axis Oz orthogonal to the plane of the disc. 
The magnetic torques are of two types: (i) A warping torque 
N m which acts in a small interaction region r ln < r < n n t, 
where some of the stellar field lines are linked to the disc 
in a quasi-cyclic fashion (involving field inflation and re- 
connection). These field lines are twisted by the differen- 
tial rotation between the star and the disc, generating a 



toroidal field AB<* 



TC-B^" 1 from the quasi-st atic verti- 
cal field Bf threading the disc, where £ ~ 1 (Aly 1985 



,(<0 



Lovelace et al.|1995 1 and the upper/lower sign refers to the 



value above/below the disc plane. Since the toroidal field 
from the stellar dipole fii is the same on both sides of 

the disc plane, the net toroidal field = + ABj, dif- 
fers above and below the disc plane, giving rise to a vertical 
force on the disc. While the mean force (averaging over the 
azimuthal direction) is zero, the uneven distribution of the 
force induces a net warping torque which tends to push the 
orientation of the disc angular momentum I away from the 
stellar spin axis (1> S (see Paper I for a simple model for this 
effect, involving a metal plane in an external magnetic field), 
(ii) A precessional torque N p which arises from the screen- 
ing of the azimuthal electric current induced in the highly 
conducting disc. This results in a difference in the radial 
component of the net magnetic field above and below the 
disc plane and therefore in a vertical force on the disc. The 
resulting torque tends to cause I to precess around £j s . In 
Paper I, we parametrized the two magnetic torques (per unit 
area) on the disc as 



N w = -(Er fi)cos/?r ro Z x (w s x I) 
N p = (Er 2 fi) cos /? f2 p w s X Z, 



(5) 
(6) 



where S(r) is the surface density, Q(r) the rotation rate of 
the disc, and j3{r) is the disc tilt angle (the angle between 
Z(r) and the spin axis a> s ). The warping rate and precession 
angular frequency at radius r are given by 
.2 



T w {r) = 



fi P (r) 



7r 2 r 7 n(r)E(r)D(r) 



-F(6>* 



(7) 

(8) 



where (9* is the angle between the magnetic dipole axis and 
the spin axis, and the dimensionless function D(r) is given 
by 



Dl 



1, ^2H(r)/r in 



(9) 



with H(r) the half-thickness of the disc. The function -F(#*) 
depends on the dielectric properties of the disc. We can write 



F(6y = 2/ cos 2 0* - sin 2 I 



(10) 



netic); see Lai & Zhang I 2008 I 



If the stellar vertical field is entirely screened out by the disc, 

3 In the notation of Long et al., r? = and = 1 corre- 

sponds to the solution for spherical accretion 
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the parameter / = 1; if only the time- varying component of 
that field is screened out, we get / = 0. In reality, / lies 
between and 1. 

The magnetic torque formulae given above contain un- 
certain parameters (e.g., f, which parametrizes the amount 
of azimuthal twist of the magnetic field threading the disc); 
this is inevitable given the complicated nature of magnetic 
field - disc interactions. Also, while the expression for the 
warping torque [eq. ||5|] is formally valid for large disc warps, 
the expression for the precessional torque was derived under 
the assumption that the disc is locally flat [eq. S is strictly 
valid only for a completely flat disc (Aly 1980)]; when this 
assumption breaks down (i.e., when \dl/d\nr\ is large), we 
expect a similar torque expression to hold, but with modi- 
fied numerical factors (e.g. the function D(r) in eq. (JsJ) will 
be different). In the application discussed in the following 
sections, we find that the condition \dl/d\nr\ <J 1 is always 
satisfied. Thus we believe that our simple formulae capture 
the qualitative behavior of accretion discs subject to mag- 
netic torques. 

It is also worth noting that the expressions |5]|6| for the 
torques only correspond to the zero-frequency component 
of the magnetic forces acting on the disc. The time varying 
components of these forces can also have significant effects. 



In particular, Lai & Zhang (2008) discussed how the com- 



ponents of the magnetic forces varying at the stellar spin 
frequency and at twice that frequency can excite bending 
waves in discs, while Terquem & Papaloizou (2000) showed 



that if the star has a dipole field misaligned with its rotation 
axis, magnetic effects create a steady-state warp in a frame 
corotating with the star. However, these "dynamical waves" 
average to zero over the stellar rotation period and do not 
affect the secular evolution of the stellar spin. In this paper, 
we concern ourselves only with long-term effects, effectively 
studying a disc profile averaged over multiple stellar rota- 
tions. 



2.2 Spin Evolution of the Star 

The effects of the magnetic torques on the evolution of the 
star - disc system are twofold. First, they will cause the 
orientation of the disc l(r) to deviate from a flat disc pro- 
file Z(r) = J ou t = J(r ut), set at the outer disc radius r out . 
These deviations will be studied in details for different disc 
parameters in Sections 3-5. Second, the back-reaction of the 
torques will change the orientation of the stellar spin axis on 
a longer timescale. The secular evolution of the stellar spin 
under the combined effects of matter accretion and star - 
disc interactions is explored in Paper I in the case of flat 
discs. Here we generalize the basic formulae derived in Pa- 
per I to warped discs. 

In general, the spin angular momentum of the star, 
J s u> s , evolves according to the equation 

~ (Js&s) = Af = Afi + Af s + Af w + Af p . (11) 
at 

Here Afi represents the torque component that is aligned 
with the inner disc axis Z(n n ) = Zi n . We parametrize Afi by 



Afi = AJlf(GM t r m ) 1/2 lm = Wo h 



(12) 



Equation ( |12[ ) includes not only the accretion torque car- 
ried by the accreting gas onto the star, M acc (GMri n ) 1 ^ 2 Z 



(where Af acc may be smaller than M, the disc accretion 
rate) , but also the magnetic braking torque associated with 
the disc - star linkage, as well as any angular momentum 
carried away by the wind from the magnetosphere bound- 
ary ( Shu et al.|1994 Romanova et al.|2009 1. All these effects 
are parametrized by the parameter A ^ 1. In particular, if 
a wind carries away most of the angular momentum of the 
inner disc, we may get A< 1. 

The term Af 3 = — \Af s \u> s represents a spindown torque 
carried by a wind/jet from the open field lines region of the 
star (e.g. Matt & Pudritz 2005). The terms Af w and Af p 
represent the back-reactions of the warping and precessional 
torques: 



(13) 



Since both N u , and N p decrease rapidly with radius (as r ), 
the integral can be carried out approximately, giving 

Af p +Af w ~ Afa \fipOJ a x Z in + n w i in X (iJ s X Z in )] , (14) 

with 

n K = -~ Lr-F(6g cos An , (15) 



3 TIT? 7 / 2 

C[l - (r in /r illt ) 3 
6r/ 7 / 2 



COS 9* COS An, 



(16) 



where cos An = lis • Zin- Note that both Afon p and Afon w 
are of order jj, 2 /rf n . 

For a fixed outer disc orientation l ou t, the inclination 
angle of the stellar spin relative to the outer disc, /3* = Aut, 
evolves according to the equation 

Js ~ COS A = Af ■ lout — cos A (Af ■ &s) 

« Afo [A (Z in ■ Zout - COS A COS An) 

+ n w (cos A — cos An Zin ■ Zout — cos A sin 2 An) 

+ 7lpU> s ■ (Z in X Zout)] ■ (17) 

Note that this does not depend on the specific form of Af 3 . 



For flat discs, equation (171 reduces to (Paper I) 
d 



dt 



cos A 



Afo 

Js 



sin 2 A (A - C cos 2 a) , (18) 



with 



c = 



C[l- (nn/nnt) 3 ]cos 2 0* 
6rf/2 



(19) 



In the flat-disc approximation, the star - disc systems can 
thus be divided in two classes with very different long-term 
spin evolution (see Fig. [l]). If C, < A, cos A always increases 
in time and the system will be driven towards the aligned 
state (A = 0). On the other hand, if f > A, there are two 
"equilibrium" misalignment angles (defined by dfi^/dt = 0): 



cos At — i 



(20) 



The smaller angle Af corresponds to a stable equilibrium, 
while A- is unstable. Thus, the final state of the systems 
depends on the initial misalignment angle A(^ — 0). If 
A(i = 0) < A-! t ne system will be driven towards a moder- 
ate misalignment Af < 90°; otherwise it will evolve towards 
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Figure 1. Time derivative of the inclination angle of the stellar 
spin (in degrees) relative to the outer disc for a flat (unwarped) 
disc. The upper panel is for £/A = 0.5 < 1, in which case the spin 
evolves towards alignment. The lower panel is for £/A = v2 > 1, 
in which case the spin either evolves toward /3+ / or toward 
jS* = 180°, depending on the initial value of ft*. The quantity 
tspin is defined in equation |3| . The arrows show the direction of 
the evolution of /?* in different regions of the parameter space. 



a completely anti-aligned configuration (ft* = 180°). From 
these results, we can see that, according to the flat-disc ap- 
proximation, if A <C 1 a misaligned configuration is strongly 
favored. 

The probability distribution of the different cases for as- 
trophysical systems will thus depend on the unknown value 
of the parameters of our model, as well as on ft*(t = 0) — 
which depends on the formation history of the star - disc 
system and is quite uncertain ( Bate et al.|2010' l. For exam- 
ple, for an isotropic distribution of Z ou t on the unit sphere 
and f > A, a fraction 0.5(1 — y^) °f * ne systems would be 
anti-aligned, while the rest would tend towards a misalign- 
ment ft+. The real distribution of disc inclination is certainly 
more complex, as A and £ will vary from system to system, 
and the distribution of the initial misalignment ft* (t — 0) is 
probably not isotropic. Additionally, the orientation of the 
outer disc might vary in time. For more details on the dis- 
tribution of final inclination angles ft*, see Paper I (Sec. 5), 
as well as Section [7] of this paper, which discusses a process 
to reach anti-alignment starting from ft* (t = 0) < ft- . 

In general, the magnetic torques i nduc e disc warping 
so that I depends on r, and equation (171 must be used 
to determine the long-term spin evolution. Since the disc 
warp evolution timescale is much shorter than the stellar 
spin evolution timescale, the steady-state warp profile l(r) 
must be solved before equation (17 1 can be applied. In the 



following sections, we will show that for most (but not all) 
realistic choices of the free parameters in our model, using 



equation (181 instead of equation (17 1 does not significantly 



change our qualitative description of the long term behavior 
of the system. We will measure deviations from the flat-disc 



approximation through the parameter £ defined by 



| C0S ^ = H| C0S ^ 



(21) 



where the left-hand side is computed using the first line of 



equation ( 17 1 



3 DESCRIPTION OF WARPED DISCS: 
THEORY 

As noted in Section 1, the evolution of the coupled star - disc 
system occurs over two different timescales. The first, tdisc, 
characterizes the evolution of the disc under the magnetic 
torques and internal stresses towards a warped steady-state 
configuration, assuming that the spin of the star u> s is fixed. 
The second, t Bp in, determines the evolution of u> a due to the 
combined effects of mass accretion and magnetic torques. 
Since we expect tdisc <C i sp in, if the orientation of the outer 
disc is fixed we can consider that, at all times, the disc is in 
a steady-state l C q(r; <2> s ). The evolution of the system is then 
described by a sequence of steady-state profiles l cq (r,u> s (t)) 



where w s (t) evolves according to equation (17 1 applied to 



Z(r) = l eq (r; ui s (t)). As discussed before, the disc itself will 
always show variations on shorter timescales (of the order of 
the stellar rotation period), which do not affect the secular 
evolution of the stellar spin and are averaged over in our 
description of the system. 

Here we describe our method to calculate the evolution 
and steady-state of warped discs. 

Systematic theoretical study on warped discs began 
with the work of |Papaloizou fc Pringfe] ( |1983[ ) and |Pa-| 
paloizou & Lin (19951, who showed that there are two dy- 



namical regimes for warp propagation in linear theory (for 
sufficiently small warps). For high viscosity Keplerian discs 
with a J> 8 = H/r (where a is the Shakura-Sunyaev pa- 



rameter so that the viscosity is 



otH SI), the warp 



satisfies a diffusion-type equation with diffusion coefficient 
U2 = v/(2a 2 ). For low- viscosity discs, on the other hand, 
the warp satisfies a wave-like equation and propagates with 
speed QH/2. In the diffusive regime, the linear theory of 
iPapaloizou & Pringle ( 1983[) was generalized to large incli- 



nation angles by Pringle| ( |1992[ ) in the limit of small local 
variations of the disc inclination. A fully nonlinear theory 
was derived by Ogilvie (19991, with prescriptions for arbi- 



trary variations of the inclination. The basic features of the 
theory were recently confirmed by the numerical simulations 
of Lodato & Price (20101. For low- viscosity Keplerian discs 
(a H/r), the linearized equations for long wavelength 



bending waves were derived by Lubow and Ogilvie ( 2000 \ 



and Lubow et al. (2002), and a theory for non- linear bend- 



ing waves was developed by Ogilvie ( 2006 \ 



For protostellar discs, recent work by Terquem ( 2008 1 



suggests that far away from the star the disc could have a 
very small viscosity parameter (a ~ 10~ 2 — 10~ 4 ), and would 
thus be described by the formalism of |Lubow and O gilvie 
(2000). However, close to the star (around a few stellar radii) 
where magnetic effects are most important and the disc warp 
can develop, the value of the effective viscosity is unknown. 
Thus in this paper, we will study both high- viscosity discs 
and low- viscosity discs. 
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3.1 High- Viscosity Discs 

3.1.1 Evolution Equations 

For viscous discs satisfying a <; H/r, we start from the equa- 
tions derived by |Og ilvic ( 1999]) . The main evolution equa- 
tions for the disc are the conservation of mass 



<9E 1 d 
dt r dr 

and angular momentum 



{tYVr) = 



(22) 



dt 



r dr 



r dr 



+ ~ I Qilr'tf^ + Q 3 Ir 3 n 2 l X ^ ) - N„, . 
r or \ or Or 



(23) 



where Vr is the average radial velocity of the fluid at a given 
radius. The coefficients Qi.2,3 characterize the magnitude of 
the various viscous interactions, while 



I = 



2tt 



dcf> 



2 j 
pz az 



(24) 



depends on the vertical density profile of the disc. The term 
N m = N u , + N p is the external magnetic torque per unit 
area. 

In general, the viscous coefficients Qi,2,3 are functions 
of the viscosity parameter a, the warp amplitude ip 2 = 
\dl/dlnr\ 2 , and the disc rotation law fl. Their values can be 
obtain ed through nu merical integration of a set of coupled 
ODEs ( |Ogilvie|l999| . In the limit ip 2 -> 0, the viscous coef- 



ficients are given by equations [141-143] of Ogilvie (1999) 



2 16a 



Qi = 
Q 2 = 
Qz = 

For ip 2 = and Q3 = 0, this is equivalent to the formal- 
ism of Pringle (1992): the viscosities v\ — v and v 2 used by 
Pringle (19921, which correspond respectively to the shear 



1 

4a 
3 

8 4 



-o(^ 2 ). 



(25) 
(26) 
(27) 



viscosity usually associated with flat discs and the viscous 
torque working against the warping of the disc, are propor- 
tional to Qi and Q 2 . The additional term Q3 was discov- 



ered by Ogilvie (19991, and contributes to the precession of 
a warped disc. For the disc configurations considered in this 
paper, the effects of finite ip 2 are small — hence, our numer- 
ical results will be computed in the limit %p 2 <C 1. However, 
we do consider the effects of non-zero Q3. 



To obtain numerical solutions to equations 1 22 1-( 23 1, 



it is convenient to switch to the logarithmic coordinate 
p — ]n(r/r- m ). We then define the logarithmic derivative 
' = d/dlnr = d/dp and the warp amplitude ip 2 = \l \ 2 . 
From equations ( |22)23[ ), we can derive the radial velocity 



V R 



Er(r 2 fi)' 



[(/Qir 2 n 2 )' 



IQ 2 r 2 n 2 ip 2 ] 



Using equation (28 1 in ( 22 l-( 23 1 then yields 



,9E 
dt 



(iQ ir 2 n 2 y - iQ 2 r 2 n 2 ip 2 



{r 2 D.)' 



(28) 



(29) 



and 



' dt ( Er2 ™) + ijkry (W 1 ^ 2 )' ^ JQ*-W) ' 

-Ir 2 tl 2 (qJ + Q 2 l' + Q 3 l x I')] ' = r 2 N m . (30) 



3.1.2 Disc Model 

For our numerical calculations, we co nsider Keplerian discs. 
If we compare the projection of (23 1 along I for I — with 
the standard flat-disc equation 

3 cKT 



di (Er n)+ rd~r 



EVBr J n - i/Er 



we see that 



Using Qi[V> 2 
have 



(32) 



Qi[V> 2 = 0]lr 2 i} 2 
= 0] = -3a/2 and Q = VGMr~ 3 , we then 



dr 
dr 



= 0, 



(31) 



E/r 



(33) 



We also rescale the time coordinate by the viscous time eval- 
uated on the inner edge of the disc: r = t/t v i s (ri n ), where 

„2 



tv 



r 

V2 



and 



;/ 2 



2Q 2 H 2 Q ~ ^-H 2 Vt 
2a 



(34) 



(35) 



is the viscosity associated with the vertical shear in the disc. 
By projecting the evolution equation of the disc angular 
momentum onto directions parallel and orthogonal to I, we 
find 



d 

d-r° 

±1 

dr 



-3/2 

Q 2 

-3/2 Ql 

oQ 2 



[s' - (S' ■ J)z] + 

C cos 2 ( 



(36) 
(37) 



p3j? 7/2 \irD(p) " 



-I x (6j s x Z) 



where the new variables u, S and p are defined by 
Er 

(rEfl at )|r=r in 



S = 



P = 



/ o Q 2 2 \ ; Q 2 ~J Q3 
a — — — -pr-cryj ] I 



Qi 



2Qi 



-al 



2Q : 



r 

Tin 



and Efl at is the surface density of a flat disc 

M M 
Z-kv 3naH 2 Q 



Efla 



(38) 

al x l' (39) 
(40) 

(41) 



Equations ( 36 1 and ( 37 1 form our model for the evolution of 



viscous discs interacting with a magnetic star. Note that as I 
is a unit vector, it only corresponds to two degrees of freedom 



in the system. Accordingly, equation (37 1 guarantees that 



dl/dr is orthogonal to I. In practice, to avoid introducing 
a preferred direction in the system (as we want to allow 
arbitrary inclination angles for the disc), we evolve all 3 
components of I, but normalize I at each timestep to the 
accumulation of numerical errors. 
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3.1.3 Steady-State Equations 



3.2 Low- Viscosity discs 



From equations | |36[ ), ( |37[ ) and ( |39[ ), it is fairly easy to derive 
the equations denning the steady-state configuratio n of the 
disc. If we set 81/ dr — and dajdr = = S' ■ I in (371, we 
obtain 



S' = 



(w s .Z) /C cos 2 6», ; j. F(6>*) „ ; 



p 3 7] 7 / 2 



nD(p) 



Equation (39 1 projected onto I gives 



and projected in the plane orthogonal to I gives 
2Qi 



5i [(S.!)i- 
>.o- L 



(42) 
(43) 

(44) 



For Q3 = 0, we thus have a set of first order differential equa- 
tions of the form U' = F(U). Given appropriate boundary 
conditions at ri n , it can easily be solved by numerical inte- 
gration. For Qz 7^ OjWe can still perform numerical integra- 

I "/ 
tion if we consider ( 44 1 as an implicit equation for I which 



has to be solved at each step of the integration algorithm. 

In practice however, the boundary conditions are im- 
posed partly at the inner edge rj n and partly at the outer 
edge r out . Indeed, we consider the orientation of the outer 
disc to be fixed 



J(?*out) — Jout 
and the mass accretion rate to be constant 
M = -2nrV R Y, 



(45) 



(46) 



(the sign is chosen so that M > for Vr < 0). We also 
impose a zero-torque boundary condition at the inner edge 

l{r in ) = (47) 

and set the surface density there to 

E(r in ) = a in E flat (r in ) (48) 

for some freely specifiable scalar <Ji n . Combining (|46l) with 



the zero-torque boundary condition and equations ( 28 1 and 
; boundar 

o-'[n B ] = 



(41 1, we obtain a simple boundary condition on a 1 at n 

1 



2' 



while ( 48 I gives the value of a [r- ln 



(49) 



(50) 



We thus have 4 boundary conditions at n n (on I , a and 
a') and 2 at r out (on I). To solve the system numeri- 
cally we use a shooting method starting at r- m . Writing 
I — (cos j3 cos 7, cos P sin 7, sin /3) and 6j s = (0, 0, 1) , we use a 
2-D Newton-Raphson method to solve for the values of /3[n n ] 
and 7[nn] leading to a solution satisfying l(r out ) = l ou t- The 
system of first-order ODEs which has to be solved at each 
iteration of the Newton-Raphson algorithm is treated using 
the 5th order StepperDopr5 method of |Press et al.| ( |2007[ ), 
and the integration is performed under the constraint \l\ = 1. 



3.2.1 Evolution Equations 

For discs with a viscosity parameter small compared to the 
thickness (a <> 5 = H/r), we can no longer use the evolution 
equations of Ogilvic (19991. In this case, disc warps propa- 
gate as bending waves. In the linear regime, the warp evo- 



lution equations were derived by Lubow and Ogilvie ( 2000 1 



1 dG 



dG 

dt 



dt r dr 

Q 2 -n 2 r 

^20, 



IxG- aQG - 



Er 3 c 2 fi dl 
4 dr' 



(51) 
(52) 



where c s = HQ Z is the disc sound speed, fl r and Q z are 
the radial epicyclic frequency and the vertical oscillation fre- 
quency associated with circular orbits at a given radius from 
the star, G is the internal torque of the disc, and E = Efl at is 
the surface density. These equations are only valid for a <J S, 
\nl - tt' 2 \ < SO. 2 and |!T2 2 - Q 2 \ <8tt 2 . In the following, we 
shall use Q r — f2 z = Q, although we verified that small devi- 
ations from these equalities do not significantly modify our 
results. 



Equations (51 1- ( 52 1 admit wave solutions. We define a 
Cartesian coordinate system so that l z ~ 1 and \lx,y\ <C 1, 
and the internal torque G acts in the a;y-plane. Consider a 



local (WKB) wave with l x>v ,G oc e % 



in a Keplerian 



disc with N m = 0. For u f2, the dispersion relation of the 



wave is, neglecting the damping term aQG in ( 52 1 

oj c s H£l 

k~ ~2~ ~2~' 
with the eigenmodes satisfying 



r 3 c s !]E 



— = F 



(53) 



(54) 



Mr 2 Q 

The + mode and — mode correspond to the outgoing and 
ingoing bending waves, respectively. 

A generic warp perturbation will not behave as pure 
eigenmodes. For numerical evolutions, it is convenient to 
define the variables 



±x,y 



— lx 



G, 



(55) 



r Mrf n fi[r in ] 

Then, the evolution equations for the disc can be written as 



d_ 



v± 



2p 3 / 2 



T V± + (V- 



4 6 



(56) 



cos_P_ 3a5 



F(6+ 



-W, X I 



C COS 2 0* ; 



I X (o> s X 1} 



ttD(p) 

Here the dimensionless time r = t8fl(r- m ) and length p = 
r/ri n are chosen so that the sound speed at the inner edge 
of the disc is c s (ri n ) = HQ. Z = 1. For the computation of 
the magnetic torque, we use the following approximations, 
accurate to first order in l x . y : 

lj s x I = — cos /3l v e x + (sin/3 + cos /3l x )e y (57) 
I x (a> s X I) = — (sin/? + cos /3l x )e x — cos /3l v e y . (58) 

The boundary conditions are particularly simple to im- 
plement for this choice of variables. At the outer edge of the 
disc, we require the ingoing mode to vanish 



V_(r out ) =0. 



(59) 
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At the inner edge, we impose the zero-torque boundary con- 
dition I — 0, G = 0, which can be written in terms of our 
evolution variables as 



V-(rfa) = V+(r in ), 
V'_(r in ) = -V' + (r in 



(60) 
(61) 



In terms of the propagation of bending waves, this corre- 
sponds to the requirement that the waves be reflected at 
the inner edge of the disc. 

3.2.2 Steady-State Warp 

The steady-state profile of low-viscosity discs can be ob- 



tained by numerical integration of equation ( 56 1 or equa- 



tions ( 51 H 52 1, by setting d/dr — 0. In practice however, the 
steady-state profile of a low-viscosity disc is nearly always 
very well approximated by a flat disc profile. The amount 
of disc warping can then be e valuated analytically. Noting 



3„2 



that T,r J c 
give 



r 3/2 



, equations (51 1-152 1 can be combined to 



dr 



, 3/2 —i 

dr 



4arN m 

' (r 3 £c§) in 



(62) 



Since N m is falling rapidly with r (N m ~ r 5 ), and 81/ dr — 
at r = ?-i n , we integrate the above equation from r- m to r: 



8 



4a (r 2 N m ) - (r 2 N„ 



dr 3 p'^ 2 (r 3 T,c 2 ) in 
Integrating from r out to n n , we then obtain 



^in ^out 



A6aN f , 



(63) 



(64) 



Using Eqs. we have 

|L - Lut| ^ j [tv is (|r u ,| + |n p |)sin(2#)] , (65) 

where t v i s = rf n /v2 is the viscous timescale for the warp, 
with V2 — c s H/ (2a). Thus, the distortion of the disk can be 
seen as arising from the warping and precessional torques 
acting over the disc during a time of order the viscous time 
scale at the inner disc edge (where the magnetic torques are 
the strongest). Projecting Eq. ( |64[ | in the direction of the 
stellar spin axis ui s and using Eqs. j5j-j6j, we have 



8 



cos An — cos /3 out = — - (t vis T w cos p sin P) . 



Since 



(^visTu; ; 



3a 2 C 



(66) 



(67) 



2v 7/2 — 

we see that as long as a 2 ( <C r) 7 ^ 2 , a condition satisfied for 
most parameters, the warp across the whole disc is small: 

6a 2 Csin(2/3) cos 2 0+ 



| An -P ut\ 



7r? 7 /2 



(68) 



For example, with rj <; 0.5, we find that for all discs |An — 
Amt| -C 1 if a 0.15. This is almost certainly true for discs 
in which bending waves can propagate. 

It is important to note that, although the approxi- 
mate analytical expression of the global disc distortion de- 
rived above is based on low- viscosity discs, our result for 
| An — Aut| is also valid for higher- viscosity discs. Indeed, in 



the linear regime and for Keplerian discs, the steady-state 
equations are identical regardless of the viscosity regime con- 
sidered. 



STEADY-STATE PROFILE OF WARPED 
DISCS AND BACK-REACTION ON 
STELLAR SPINS 



Using the numerical scheme presented in Section [3.1.3[ we 
can now determine the time-averaged steady-state profile 
of the disc under the influence of the torques exerted by a 
magnetic star. The characteristics of the warped disc will of 
course vary with the choice of the free parameters included 
in our theoretical model. We begin our study by showing 
results for two standard discs, chosen so that they belong 
to the two classes of long term stellar spin evolution pre- 
dicted in Section [2.2| when the accretion parameter defined 
in equation ( 12 1 is A ~ 0.5 (a typical value in the allowed 



range ^ A 1). We then vary the disc parameters, and dis- 
cuss their influence on the disc profile, and on the spin evo- 
lution. Finally, we check that, as predicted in section ["3. 2. 2| 
low- viscosity discs, which follow the different evolution equa- 
tions described in section 3.2.1 (valid for a <J S = H/R) have 



only negligible steady-state warps and are for all practical 
purposes well described by the flat-disc approximation. 

Our base models are discs with viscosity a = 0.15 and 
thickness 8 = 0.1. We fix the surface density at the inner 
boundary by setting a ln = 1.0 so that it is equal to the 
surface density of a flat disc ( |41[ ), choose the inclination an- 
gle of the outer disc A>ut = P* = 10° and the magnetic 
inclination angle 8+ = 30° with respect to the spin <2} s . 
The star is assumed to have mass M* = M and radius 
Ri, — 2Rq The strength of the magnetic field is chosen so 
that n n = 2.5-R* (Corresponding to ~ lkG for typical 
parameters, see Eq. [2]), and the action of the torque N m is 
limited to the region n n ^ r ^ ri Ilt = 1.5n n . The accretion 
rate is M = 10 -8 Mo/yr, and we put the outer disc bound- 
ary at r out = 10 4 ri n (corresponding to r out w 250AU, a size 
typical of the observed protoplanetary discs). This disc has 
small values of ip 2 everywhere, and accordingly we neglect 
the nonlinear terms in Qi (but we keep Qz = 3/8). The other 
parameters are chosen to be Q — 1, / = 0, and either rj = 1 
(so that the long-term evolution of the system aligns the 
spin axis lj s with the disc axis) or r\ = 0.5 (for which the 
flat-disc approximation predicts a long term misalignment 
toward Af ~ 45° if the initial disc has A ^ 135°). 

We should note that these parameters are purposefully 
chosen to test the limits of the flat disc approximation. Our 
choices M*, i?*, B t , S and M© are relatively standard values 



for protoplanetary discs around T-Tauri stars [see Bouvier 
et al. (2007b I and references therein], while there are no 



particular reasons to prefer any specific orientation of the 
magnetic dipole 6+. But a = 0.15 is larger that recent esti- 



mates of the viscosity in the outer parts of the disc ( Terquem 



2008L and probably on the high end of what can be expected 
in the inner disc. However, we have shown that smaller val- 
ues of a lead to smaller amplitudes of the steady-state warp 
(the warp amplitude is proportional to a 2 ). Thus, the flat 
disc approximation is more likely to be satisfied at low vis- 
cosities. 

In order to analyze the radial variations of the disc warp 
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profile, we define the tilt j3[r] and the twist y[r] by 

l[r] — (sin /3[r] cos 7[r], sin/3[r] sin7[r], cos /3[r]), (69) 

with the convention that 7[r out ] = 0. Some parameters of 
the system can be varied without modifying the dimension- 
less solution for the profile of the surface density cr(p) and 
the orientation of the disc l(p): modifications of M, M*, i?* 
or r- m (at constant 77, r ou t/r- m and ri nt /n n ) will influence the 
values of the timescales t v is and t apm , but not f3[p] or 7[p]. 
Thus, the steady-state profile can be solved while keeping 
these parameters fixed without any loss of generality. The 
disc profile in physical units [S(r), l(r)] can easily be re- 
trieved from the dimensionless solution [cr(p), l(p)]. Addi- 
tionally, the four parameters (77, 6**, /, f) correspond to only 
two degrees of freedom in the model, through the quantities 



equations (64 1 and (681. In particular, formula (68 1 applied 



ci 



C2 = 



C cos 2 0* 
rj 3 - 5 

F(0*) = 

„3.5 



2/ cos 2 



(70) 
(71) 



We will thus limit ourselves to variations of £ and /. Varying 
the thickness 8 of the disc has very similar effects: it changes 
the value of the function D(r) at small radii, effectively mod- 
ifying the value of C2 close to n n . As the magnetic torques 
mostly affect the region close to the inner edge of the disc, 
the influence of 8 is similar to that of F(8+). Finally, we 
are also free to modify the boundary conditions used, and 
in particular the choices of r out and o"i n . Varying r 

out seems 

to have only negligible effects, as long as r out /ri n is large 
enough for a steady-state solution to exist. Decreasing o"i n , 
on the other hand, leads to more significant changes in the 
warp profile. A small cri n favors warping disc, so that a de- 
crease of o"i n has an effect similar to increasing both c\ and 

C2- 

Thus, the effects of varying various parameters of the 
system can be examined with our standard discs, by varying 
only two parameters, £ (or ci) and / (or C2). In section 4.1 
we present our results for our two standard discs, which are 
similar except for the value of the parameter 77 (changing 
77 correspond to a rescaling of both the warping and the 



precessional torque). Then, in section 4.2 we study varia- 



tions of the warping torque alone, by modifying the value of 
the parameter ( characterizing the strength of the toroidal 
field in the disc. The influence of the precessional torque is 
studied in more details in section |4~3| through variations of 
the parameter / (related to the ability of the time-varying 
component of the vertical magnetic field to penetrate the 



disc). Finally, in section 4.4 we comment on the influence of 



the parameter Q3, which was usually neglected in previous 
studies of warped discs. 



4.1 Standard disc results 

The profile for the tilt and twist angles of our standard con- 
figurations (Fig. [2| show a relatively weak warping of the 
disc. For the disc with weaker magnetic interactions (77 = 1), 
the difference in tilt between the inner and outer edges is 
about 0.17° and the twist over the whole disc is 1.2°, while 
for stronger interactions (77 = 0.5) the disc is tilted by 1.8° 
and twisted over 16°. These warps are comparable in magni- 
tude to what we could have predicted using the approximate 



to these two choices of parameters predicts tilts of 0.28° and 
3.2°, respectively, with most of the difference between the 
approximate formula and the numerical results due to the 
cutoff applied to the magnetic torques at r = n nt , neglected 



in the derivation of ( 68 1 



Note that we choose to vary the parameter 77 defined in 
equation Q as it conveniently modifies the effective strength 
of both magnetic torques in our model. In practice, 77 is de- 
termined by the geometry of the accretion flow, while un- 
known physical parameters such as the dipole strength p, its 
orientation (9*, the surface density at the inner edge of the 
disc am or the magnetic twist parameter £ will vary from 
system to system. 

Given the small warp, the evolution of the misalign- 
ment angl e 3± between u> s and l ou t is well approximated by 
equation (17 1 with ii„ = lout' if we compute J\f from the 
steady-state profile Z(r), we find that the parameter £ in 
equation (21 1, which parametrizes deviations from the flat 
disc approximation (£ = 1 for a flat disc) is £ = 0.997 for 
rj — 1 and £ = 0.86 for 77 = 0.5, if we set the accretion 
parameter A to (we choose A = when computing £ in or- 
der to measure directly differences in the effect of the back- 
reaction magnetic torques between the flat-disc model and 
the warped disc steady-state, disentangled from the effect of 
angular momentum accretion). 

However, even a small disc warp can significantly change 
the critical angles /3± for which dj3*/dt = 0. By varying 
/3* = /3(r ou t), we can determine the values of /3± numer- 
ically. For 77 = 0.5 and A = 0.5, we find /3+ = 32° and 
/3_ = 148° which are quite different from the prediction of 
the flat-disc approximation (/3+ = 45°, /3_ = 135°). This is 
due mainly due to the effect of the twist of the disc. In the 
flat-disc approximation, 7 = and the back-reaction due to 
the precession torque has no effect on the evolution of the 
stellar spin. But as long as 7i n < 7T, that back-reaction will 
tend to align the stellar spin and the disc orbital angular mo- 
mentum, and this effect can be large enough to significantly 
shift the value of /3± (see also subsection |4.3| . The same 
effect can also modify the qualitative behavior of systems 
for which the predicted misalignment angle /3+ is close to 0, 
in such a way that the stable configuration at /3+ no longer 
exist. The orbital angular momentum of the disc would then 
be expected to align with the direction of the stellar spin. 



4.2 Large warping torques 

Realistic discs are expected to have parameters £ ~ 1 (char- 
acterizing the azimuthal magnetic twist) and 77 ~ 0.5 (char- 
acterizing the inner disc radius; see Eq. [I]), but the exact 
values of those parameters are unknown (see Section [2J . By 
increasing £ we see that the approximation i(r; n ) = Z(r ou t) 
can break down for high-viscosity discs. In Fig. [3] we show 
the variation of the steady-state disc profile when £ is varied 
between 1 and 5 for 77 = 0.5 and our standard disc parame- 
ters. Clearly, there can be large differences between the ori- 
entation of the disc at its inner and outer edges when £ <; 1. 
In Fig. [4] we show the value of £ [Eq. [2l] for various choices 
of £. At low £ <y 4 and for the choice of accretion parameter 
A = 0, the flat-disc approximation predicts the magnitude of 
the back-reaction torques acting on the star within a factor 
of 2. Deviations at low £ w 0.5 are due to the relatively large 
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Figure 2. Upper panel: Disc tilt angle ft for our standard disc, 
with r\ = 1 and »; = 0.5. Lower panel: Twist angle 7 for the same 
parameters. 




log(r/r. n ) 

Figure 3. Upper panel: Disc tilt angle /3 for different choices of £ 
values, all with r\ = 0.5. Lower panel: Twist angle 7 for the same 
disc parameters. The outer edge of the disc is at r ou t = 10 4 r; n . 



influence of the precessional torque (which is independent of 
when the warping torque becomes small. The long-term 
evolution of the stellar spin direction will remain similar to 
the flat-disc predictions, with a stable configuration at some 
misalignment angle (3+ 7^ for most values of the accre- 
tion parameter A. But for £ 4.5, the twist is so large that 
the behavior is the opposite of what would be predicted by 
our approximate flat-disc formula: the back-reaction tends 
to align the disc and the spin of the star. This shows that 
at large £ we must determine for each set of parameters 
the profiles /3[r] and y[r] in order to predict the long term 
evolution of the stellar spin. As an example, we construct 
sequences of steady-state disc configurations for a fixed f, 
varying the inclination angle of the outer disc /3*. In Figs. 
5|7 we show the resulting d cos /3*/dt for £=1,3 and 5, and 
compare with the predictions of the flat-disc approximation. 
For £ = 1,3, the general behavior is similar to what the flat- 
disc approximation predicts. As seen in subsection |4.1[ the 
precessional torque will favor alignment of the stellar spin 
with the disc orbital angular momentum, so that the numer- 
ical results usually show that /?+,jv um ^ P+,Flat — at least 
as long as F(0+) is of order unity. For £ = 5, however, signif- 
icant differences become visible. At small inclination angles, 
the system will evolve towards /3* = 0, while at large incli- 
nations, the system will evolve towards /3* « 165°. In the 
intermediate region 15° ft* Ss 135°, the system will evolve 
towards /3* w 50° (for A = 0.5). Finally, for larger £ we 
are in a completely different regime: for some inclinations, 
two steady-state solutions exist. Clearly, to determine which 
of those steady-state solution is relevant requires numerical 
integration of the time evolution of the star-disc system. 

Our current understanding of the effects of magnetic 
fields close to the inner edge of the accretion disc is not 
sufficient to determine with certainty the range of realistic 
values of the parameters £ and r\. However, their favored 
values lie in a region of parameter space where the flat disc 




Figure 4. Variation of the parameter £ characterizing the devia- 
tions from the flat-disc approximation [see Eq. | |21[ l] as a function 
of £, for a sequence of discs with r\ = 0.5 and the choice A = for 
the accretion parameter [see Eq. l |f 2| ], 



approximation appears to hold relatively well (f ~ 1, 77 ~ 
0.5). The large deviations from the flat disc model observed 
at high C, are thus unlikely to be encountered in astrophysical 
systems, though they cannot be entirely ruled out. Thus, 
these results implies that the flat disc approximation is likely 
to be justified, with the caveat that it tends to overestimate 
the value of the misalignment angle /3+. 



Spin 
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Figure 5. Secular evolution rate of the spin-disc inclination angle 
0* for discs with £ = 1 . The time derivative of cos /3* is given for 
the flat-disc approximation (Flat) and for our numerical results 
for warped discs (Num), as well as for 3 different values of the 
accretion parameter A = 0,0.5, l(see equation |12[ ), The angle /3* 
will increase if d cos /3* / dt < 0. 
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Figure 6. Same as Fig. [3] except that we use £ = 3. 



4.3 Varying the precessional torque 

The results presented in previous subsections were all ob- 
tained with / = 0, thus fixing the choice of the function 
F(9 t ) characterizing the magnetically driven disc precession 
rate. Using different values of /, even at low £ it is possible 
to find discs which require numerical solutions to determine 
their warp profiles. For example, if we choose / = 1 in- 
stead of / = 0, the sign and magnitude of fl p will change. 
The twist of the disc becomes more important, so that even 
for £ = 1, jy = 0.5, there is a significant deviation from 
the behavior of the flat-disc configuration. Comparisons be- 
tween the disc profiles for different / can be found in Fig. 
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Figure 7. Same as Fig. [5] except that we use f = 5. Note that 
when the disc is nearly aligned or nearly anti-aligned, the quali- 
tative behavior of the solution is different from the predictions of 
the flat-disc approximation. 

[8] The most important feature of these profiles is that, for 
the larger values of /, we have a large twist 7(n n ). Hence, 
the precession term in equation \17\ (proportional to n p ), 
which does not contribute to the evolution of /3* in the flat- 
disc approximation, now has a significant impact. For a twist 
7 such that sin(7 — y[r out ])F(9*) > 0, the precession term 
directly contributes to the alignment of the outer disc axis 
with the stellar spin. This is always the case for disc twists 
7in| ^ 180°, as a positive F(6+) causes the inner disc to 
precess in the prograde direction, while F(0+) ^ causes 
a retrograde precession. If the precessional torque becomes 
large enough compared to the warping torque (proportional 
to n w ), the long-term evolution of the stellar spin direc- 
tion will be modified. For our standard parameters and the 
choices of 77 = 0.5 and A = 0.5, we find that discs with 
/ J; 0.5 will always lead to spin-disc alignment, contradict- 
ing the flat disc predictions (see Fig.|9|. However, it is worth 
noting that some configurations with high / still allow for 
long term misalignments: for example, for / = 0.5, increas- 
ing the strength of the azimuthal B-field to £ = 3 leads 
to a behavior very similar to what we found for / = 0, 
£ = 3 (see Fig. [6JI, while decreasing the viscosity parameter 
to a = 0.015 (and choosing 8 — 0.01) limits the twist of the 
disc, so that the flat-disc approximation remains valid. 

The above results show that the precessional torque can 
in principle cause non-negligible deviations from the flat- 
disc model. Nevertheless, for the largest part of the favored 
parameter space (small a, or large a with small precessional 
torque), the flat-disc approximation is justfied. 

4.4 Influence of Q3 

As mentioned before, most previous works on warped discs 



have been done using the formalism of Pringle ( 1992 1, which 
corresponds to Qz — in the formalism of Ogilvie (19991. 



This is a good approximation, as long as the influence of the 
small precessional torque due to Q3 7^ is negligible. For 
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Figure 9. Same as Fig. [5] except that we set A = 0.5 and choose 
/ = 0, 0.16, 0.5, 1. For / = 0.16, the disc twist is very small, and 
our numerical result matches the flat-disc approximation better 
than for / = 0. For larger values of /, the disc twist is large, and 
the disc will align with the stellar spin regardless of the initial 
value of B+. 



the system studied here a small change in the twist of the 
disc can affect whether a configuration will align over time, 
or be driven towards a stable misaligned steady-state. In 
Fig. |10[ we show the difference in the disc tilt and twist for 
our standard disc with r/ = 0.5, using both the formalism of 



Ogilvie ( 1999 \ and Pringle ( 1992 1. Differences in the warp of 



the disc of a few degrees are observed, though the warps are 
small in both cases. Because the precessional torque acting 
on the disc using Q% = is smaller, it will be less twisted. 



This leads to a behavior slightly closer to what the flat-disc 
approximation predicts. If we choose the accretion parame- 



ter A = 0.5 (equation 121, then the flat-disc approximation 
predicts a stable misaligned configuration at f3 + = 45°. For 
warp discs, we find that the misalignment angle is signif- 
icantly smaller, /3+ = 32°. The difference in /3+ between 
profiles obtained using Q3 = 3/8 and Q3 = is only 0.5°, 
which is negligible at the level of accuracy our model can 
achieve. 

Q3 can also influence the qualitative behavior of the 
steady-state solutions at high-£. For strongly warped discs, 
it is sometimes possible to have two solutions satisfying the 
steady-state equations. Choosing Q3 7^ seems to limit the 
size of the region of parameter space where this happens. 
For example, for / = 0, r\ = 0.5, 9* = 10° and £ = 5.5, 
two profiles are acceptable steady-state solutions if we chose 
Qz — 0, while for Q3 = 3/8 the same parameters lead to a 
unique solution (see Fig |ll[ ). 



4.5 Low- Viscosity Discs 

In the linear regime, the equations determining the steady- 
state profile of the disc are identical for the a ^ 8 and 
a ^ 8 cases. In the previous subsections, we have seen 
that our approximate formulae for the amplitude of the 



warp, equations (64 1 and (68 1, give relatively good results 
for a ~ 8 = 0.1. We also confirmed numerically the a 2 de- 
pendence of the warp of the disc, shown in Figure |12| For 
smaller viscosities, a ^ 8, we expect the warp to be even 
smaller, and the linear approximation more accurate. Hence, 
we can immediately deduce that the time-averaged warp of 
low-viscosity discs will be extremely small. For such discs, 
the flat-disc approximation will nearly always give accurate 
results for the secular evolution of the stellar spin. 
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Figure 13. Time evolution of the disc tilt angle profile for stan- 
dard discs with a = 0.15 and r out = 100r; n . Time is in units of 
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Figure 12. Disc tilt angle /3 for a = 0.15, 0.015, 0.0015. To check 
the a 2 dependance of /3, the deviation from a flat disc is multiplied 
by 10 2 and 10 4 for a = 0.015 and a = 0.0015 respectively. 



5 TIME EVOLUTION OF DISC WARP 
TOWARD STEADY-STATE 

Having established the steady-state of warped discs, we now 
study their time evolution starting from some generic ini- 
tial conditions, when the symmetry axis of the outer disc 
is misaligned with the stellar spin. To this end, we evolve 



equations ( 36 1-( 37 ) for high- viscosity discs and (561 for low- 



viscosity discs. Since the timescale to reach steady state is 
generally much longer than the local disc warp/precession 
time F" 1 ~ fi" 1 [see Eq. nip], an implicit evolution scheme 



is necessary. Our numerical method is detailed in the Ap- 
pendix. 



5.1 High- Viscosity Discs 

For viscous discs with a *5 5 = H /r, we expect the evolution 
of the system to occur over the timescale t v - 1B (r) — r 2 /u2 = 
2a/(8 2 Q.). Note that at the disc inner edge, i V is(?"in) = 
(3a 2 Ccos 2 6'*/2j7 3 ' 5 )r^ 1 (ri n ) [see Eq. |67|] is smaller than 
the warping timescale for typical parameters. In terms of 
the dimensionless time r = t/t v i s (nn), we expect the disc 
to reach the steady-state profile at radius r within a time 
of order r ~ (r/ri n ) 3 ^ 2 (assuming constant 8). To test this 
expectation, we evolve our standard disc model (see Section 
4) for 77 = 0.5 and different locations of the outer radius 
(''out = 100n n and r out = 1000n n ), as well as for a more 
viscous disc with a — 0.3. The disc is initialized in a flat 
configuration with I — l ou t and we observe its evolution to- 
wards the steady-state profile. In Figs. |13]15| we plot the 
disc warp profiles at times r = 10 3n//4 for n = 0, 1, ...,4 — 
by which point the viscous forces should have brought the 



t/2 



''in 



disc into its steady state up to radius r ~ 10' 

In all cases, we see that the evolution occurs approx- 
imately on the expected timescales: the local distortion of 
the disc (i.e. dl/dlnr) up to radius r does not vary much 
past the viscous timescale at that radius. The orientation of 
the disc (I), on the other hand, continues to change to ac- 
commodate the evolution of the disc at larger radii. Overall, 
the disc will reach its equilibrium profile within the viscous 



timescale t v is(f\i 



where r war p is defined as the largest 



radius at which the warp \dl/dlnr\ is significant. 

For the two simulations with the outer disc boundary at 
Tout = 100ri n , we find that r warp ~ r out and the disc reaches 
its steady-state profile within r ~ 1000. At later times, the 
evolution of the profiles becomes negligible. For the larger 
disc (r out = 1000n n ), the situation is slightly different. At 
r — 1000, the disc has reached its steady-state distortion 
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Figure 15. Same as Fig. |13| except for a = 0.3. 



up to r — 100ri n . The disc will still evolve up to r ~ 10 4 ' 5 , 
but as the warp is very small for r <; 100ri n , the changes in 
the profile are minimal. As most discs studied in this paper 



show negligible warps for r>(f0 — 10 )r ; 
steady-state to be reached within at most 



iv is (10n n ) ~ 10 4j - 



500 



Q 

015 



we expect the 



5 

OA 



yrs, 
(72) 

regardless of the outer radius of the disc. As this is much 
smaller than the evolution timescale for the spin of the star, 
we are justified to consider only the steady-state configura- 
tion of the disc when attempting to determine the long-term 
evolution of the misalignment between the stellar spin and 
the orientation of the outer disc. 
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Figure 16. Evolution of the disc tilt angle profile for discs with 
a = 0.01 and r ou t = 100ri n . The unit of time is (5n(nn)) — 1 . 



5.2 Low- Viscosity Discs 

The evolution of low-viscosity discs (a <J S) is qualita- 
tively different from high- viscosity discs. According to equa- 
tion ( |56| l , perturbations around the steady-state propagates 
as bending waves, at roughly half the local sound speed. 
Thus, we expect the disc to settle to an equilibrium within 
the propagation timescale of these waves, 



2dr 



(2 x 10 yrs) 



4 

3<5fi(r ou t) 
0.1 ( r out 



(73) 



3/2 



5 V100AU/ 

In Figures |T6l[T8l we show the evolution of the disc tilt pro- 
file j3 as the bending wave propagates across the disc, using 
our standard choice of parameters for the magnetic torques. 
We consider different discs: the first two use a = 0.01 and 
have their outer boundaries at r out = 100ri n (Fig. 16 1 and 
Tout = 1000n n (Fig. 171. The third has a higher viscosity 
a = 0.05, and r out = 1000r in (Fig. [l8|. All three simula- 



tions are started from a flat disc configuration, and show 
the same behavior: the magnetic torques perturb the in- 
ner disc, and the perturbation propagates outwards over the 
timescale t wavc . Again, this timescale is much less than the 
spin evolution timescale. 

The location of the outer disc radius apparently does 
not have a significant influence on the final state of the sys- 
tem. The viscosity, on the other hand, affects the disc warp 



amplitude as predicted by equations ( 64 \ and ( 68 1: the warp 



is proportional to a , with the amplitude |/3 ou t — An| given 
by Eq. (1681) to within a factor of two. 



6 VARIATIONS OF THE OUTER DISC 
ORIENTATION 

In the previous section, we have studied the time evolution 
of warped discs under the assumption that the orientation 
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Figure 17. Same as Fig. [16] except for r ou t = 1000n n . 
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Figure 19. Time evolution of the disc tilt angle profile ft for 
a = 0.15 and r ou t = 1000ri n , when the outer disc orientation is 
changed from /3(r out ) = 10° at t = t — At to /3(r out ) = 30° at 
t = to, with At = 10 3 t v ia(nn)- Time is in units of tvis^in)- 



6.1 High- Viscosity Discs 

We first consider a viscous disc with a = 0.15 and r out = 
lOOOYin. We choose to vary the outer disc orientation over 
At = 1000tvis(Vm) ~ tvis(100>i n ). As in the case of the evo- 
lution towards the steady-state, the evolution of the disc 
occurs on the viscous timescale t v is (see Fig. 19 1. However, 



as significant changes now take place at the outer radius, 
the new steady-state configuration will be reached in a time 
of order the viscous timescale at the outer radius t v is(fout), 
whis is larger than t v i s (»Varp) (see Section 5.1 1. Nevertheless, 
even though the steady-state is likely to be reached over a 
longer timescale than when the outer orientation is fixed, we 
still expect t v i E (r out ) to be significantly less than the evolu- 
tion time for the stellar spin t sp in- Thus, if the variation of 
the orientation of the outer disc occurs on a timescale shorter 
than tspin, the evolution of the stellar spin is well described 
by the approximation in which the disc is assumed to be in 
its steady-state configuration at all times, and adapting in- 
stantaneously to modifications of its orientation at the outer 
boundary. 



6.2 Low- Viscosity Discs 



of the outer disc is fixed. However, a protoplanetary disc is 
formed inside the star forming core of a turbulent molecular 
cloud 



[e.g. 



McKee fc Ostriker (20071], Thus in general we 



expect the outer orientation of protoplanetary discs to have 
some variations in time. In this section, we study how the 
warped disc and particularly the inner disc orientation re- 
spond when the outer disc orientation varies by some finite 
amplitude (chosen to be 20°) over a period of time short 
compared to the evolution timescale of the disc, and how 
such variations affect the secular evolution of the stellar spin 
direction. 



The same type of evolution can also be studied for low- 
viscosity discs. If we choose the viscosity parameter a — 
0.01, and change the orientation of the disc by 20° over a 
timescale At = t W avc(100ri n ), where t W ave('") is defined by 



equation ( 73 1 with r ou t replaced by r, we obtain the evolu- 



tion shown in Fig.|20| We see that a bending wave created at 
the outer boundary propagates inward, until it reaches the 
inner edge of the disc where it is reflected. The total time re- 
quired for the disc to reach a new steady-state is thus twice 
the crossing time of the bending wave, ~ 8/(3<5f2(r ut))- For 
low-viscosity discs, the condition for the steady-state ap- 
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Figure 20. Same as Fig. |19| except for a disc with a = 0.01 and 
Tout = 1000ri n , and the outer disc orientation varies by 20° over 
At = i W avc(100r in ). 



proximation to be valid when the orientation of the outer 
disc is allowed to change over time is thus 

8 



spin Z 5^ r « 4000yrs. 



(74) 



As t s 



■t, lMyr [see Eq. j3j], this condition is easily sat- 
isfied. Also note that the evolution equations of bending 
waves adopted in our analysis are based on the flat-disc ap- 
proximation. When the outer disc boundary evolves as fast 
as shown in Fig. |20| this approximation is no longer valid. 
Thus in practice, we should also require At ^> 8/[3<5S7(r out )]. 



7 APPLICATION: ANTI- ALIGNED 
EXOPLANETARY ORBITS 

Our calculations in Sections 3-5 show that for the most likely 
physical parameters that characterize a magnetic star - disc 
system, the disc warp is small. Therefore the long-term evo- 
lution of the stellar spin is generally well-described by equa- 
tion ( |18[ ), as long as the orientation of the outer disc is kept 
constant. According to (18 1, three types of spin evolution 



trend are possible, depending on the parameters of the sys- 
tem and the initial conditions (Paper I). If £ < A, the stellar 
spin and the disc axis will always align (given enough time) 
regardless of their initial relative inclination. If £ > A, mis- 
alignment between the disc and the stellar spin will develop, 
evolving towards one of the two possible final states: either 
Pi, = /3+ < 90°, or a perfectly anti-aligned configuration. 
The second configuration can only be reached if the initial 
disc has a retrograde rotation with respect to the stellar spin, 
with /3(t = 0) > 180° — /3+ = P-. In this case, to explain 
the observed expolanetary systems with retrograde orbits 
relative to the stellar spin ( |Triaud et aT7||2010| >, we have to 
require that the disc rotates in a very different direction from 
the stellar rotation axis during the time of planet formation. 
As discussed in paper I [called scenario (2) in Section 5 of 



Paper I], this is certainly possible if we consider the com- 
plex nature of star formation in molecular clouds and in star 



clusters [see also Bate et al. ( 2010 1]. 

In Paper I, we describe another potential pathway to 
create retrograde exoplanetary systems [called scenario (1)] 
starting from prograde-rotating discs. If the disc axis and 
stellar spin axis are initially nearly (but not perfectly) 
aligned, and the magnetic torques are such that the aligned 
configuration is unstable, then the misalignment angle will 
tend towards /?+, and no retrograde planets can be pro- 
duced. However, this is only true if the orientation of the 
outer disc does not vary. If instead we assume that the outer 
disc experiences a change of its orientation A/3 > /3_ — f3 + 
over a timescale At sufficiently long that this change can 
propagate to the inner disc, but short enough that the stel- 
lar spin direction does not significantly evolve over At, then 
the star - disc inclination can jump to /3 > /3_ , and continue 
to evolve towards anti-alignment. These conditions can be 
summarized as: 



A/3 > P- - 13+ = 180° - 2cos _1 



^disc ^ ^spiriT 



(75) 
(76) 



where the disc warp evolution time tdi sc ~ t v i B (r ou t) if a ^ S 
(high-viscosity disc) and t disc ~ t W av C for a JS S (low- 
viscosity disc). As we have seen in Sections 6.1-6.2, the sec- 
ond and third conditions are fairly easy to satisfy, as tdi sc is 
at most of order 10 4 yrs for a viscous disc with r out ~ 10 4 n n 
(and tdisc would be significantly shorter for a smaller outer 
disc radius), while t sp i n ~ 10 6 yrs for typical parameters [See 
Eq. q3p] - The potential to satisfy the first condition, on the 
other hand, will depend on the fraction of the disc angular 
momentum which is accreted by the star (the parameter A 
in equation 12 1 and the magnetic warp efficiency (the pa- 



rameter Q. If the star only accretes a small fraction of the 
angular momentum (A <C 1), then the angles j3± are both 
close to 90° , and small variations of the outer disc are suffi- 
cient to allow the system to jump to the retrograde state and 
eventually evolve towards the anti-aligned configuration. 



8 DISCUSSION 

The main finding of our paper is that although magnetic 
interactions between a protostar and its disc have a strong 
tendency to induce warping in the inner disc region, inter- 
nal stresses in the disc tend to suppress the warping un- 
der most circumstances. The result is that in steady-state, 
the whole protoplanetary disc approximately lies in a single 
plane, which is determined by the disc angular momentum 
at large radii (averaging out the dynamical warps which 
vary on timescales of order the stellar rotation period — 
such dynamical warps do not affect the secular evolution 
of the stellar spin). The reason for the small steady-state 
disc warp is that the effective viscosity acting to suppress 
disc warp, V2 — vi/(2a 2 ), is much larger than the viscos- 
ity (y\ = aHcs) responsible for angular momentum transfer 
within the disc ( |Papaloizou fc Pringle]|1983| |Ogilvie||1999 1. 
In fact, our anaylsis of the steady-state magnetically driven 
disc warp shows that, in the linear regime, the disc inclina- 
tion angle (relative to the stellar spin axis) varies from the 
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outer disc to inner disc by the amount [see Eqs. (661 and 
168)] 



|/5in ^out | ^ (^vis-Tu; sin 2/3) in 



a 2 Csin(2/3 in ) 



'/ 



7/2 



(77) 



where t v i B = r 2 jv2 is the viscous time and T w is the warping 
rate due to the magnetic torque. This result is valid regard- 
less of whether the warp perturbations propagate diffusively 
(for a <; H/r, high- viscosity discs) or as bending waves (for 
ct Ss H/r, low- viscosity discs). Thus, for the preferred values 
of the parameters rj ~ 0.5, £ ~ 1, we find |/?m — /3 out | -C 1 
for q <C 0.3. Moreover, our analysis of the time evolution 
of warped discs shows that, starting from a generic initial 
condition, the steady-state can be reached quickly, on a 
timescale shorter than the characteristic timescale for the 
evolution of the stellar spin orientation. 

Overall, our study of magnetically driven warped discs 
presented in this paper justifies the approximate analysis 
(based on the flat-disc approximation) of the long-term evo- 
lution of spin-disc misalignment presented in Paper I. Nev- 
ertheless, we note that even relatively small disc warps can 
modify the "equilibrium" spin - disc inclination angles /3± 
(see Fig. 1) from the flat-disc values, thereby affecting the 
"attractors" of the long-term evolution of the spin - disc 
inclination angle. If we allow for more extreme parameters 
(but still reasonable by physical considerations) for the disc 
- star systems, much larger disc warps become possible and 
qualitatively different evolutionary trends for /3 may be pro- 
duced (see Figs. [5][7] and |9|. 

Taken together, the results of this paper and paper I 
demonstrate that at the end of the first stage of the plane- 
tary system formation (see Section 1), the inclination angle 
between the stellar spin and the angular momentum axis of 
the planetary orbit may have a wide range of values, includ- 
ing alignment and ant i- alignment (see also section [Tf. Dy- 
namical processes (e.g., planet-planet scatterings and Kozai 
interactions) in the second stage, if they exist, would fur- 
ther change the spin - orbit misalignment angle. More work 
is needed to determine the relative importance of the two 
stages in shaping the properties of planetary systems. Cur- 
rently, the orbital eccentricity distribution of exoplanetary 
systems suggests that the second stage is important (e.g., Ju- 
ric & Tremaine 2008) . On the other hand, as noted in paper 
I, the 7° misalignment between the ecliptic plane of the so- 
lar system and the sun's equatorial plane may be explained 
by the magnetically driven misalignment effect studied in 
this paper. Also, the recent discovery of Kepler-9 (Holman 
et al. 2010), a planetary system with two or three planets 
that lie in the same orbital plane, seems to suggest that at 
least some planetary systems are formed in a "quiet" manner 
without violent multi-body interactions. Obviously, measur- 
ing the stellar obliquity of such "quiet" systems would be 
most valuable. 
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APPENDIX A: NUMERICAL METHOD FOR 
SOLVING WARP EVOLUTIONS 



We evolve equations ( 36)37 1 for viscous discs and ( 56 \ for 
low-viscosity discs with an implicit Crank-Nicholson evolu- 
tion algorithm inspired by the method used by Pfeiffer & 
Lai ( 2004 I to study the behavior of a disc accreting onto a 



magnetic star when the orientation of the outer disc l(r ont ) 
is aligned with the stellar spin a> s . The evolution equations 
are all of the form 

d d 2 d 

^y = A{x,y)-^y + B(x,y)—y + C(x,y), (Al) 

and are discretized at the N vertices ffio,i,...,iV-i of our nu- 
merical grid as 



d_ 

Or 



d 2 



- 

At 

7 Vi+i + Vi-i - 2y« 
Ai 2A^ ~ + 



(A2) 
(A3) 



A, 



y%+i + yi-i - 2y t 
2Ax 2 



£, Vi+1 — Vi-1 D Vi+1 — Vi-1 , . .s 
a i 77 r "i 77 I A4 ) 

4.Ax 4.Ax 



C(x,y) 



-ACi + d) 



(A5) 



where yi is the value of y at point Xi and time t+At. At each 
time step of the Crank-Nicholson algorithm, we start from 
an initial guess for yi obtained by extrapolating from the 
three previous time steps. From that guess yf , we evaluate 
A, B and C. Assuming these functions as fixed, we can then 
obtain yt by solving a tridiagonal system of equations. This 
gives us an updated guess y\ for the value of the function 
at t + At. We then repeat the operation until the step s for 
which the condition 



max \ yi 



Vi 



< ttr 



(A6) 



is satisfied for some chosen tolerance e t ri- 

The main advantage of this implicit method is that the 
time step At can be much larger than the Courant limit 
when the variable y evolves slowly in time. In practice, At 
is chosen so that the condition 



max|y, - y t \ < ecN 



(A7) 



is satisfied for e tr i <S ecn -C 1. We choose ecn ~ 10~ 4 in our 
simulations (the Crank-Nicholson algorithm is second-order 
convergent in time, and we verified both the convergence 
and the fact that we could obtain sufficient precision for 
that choice of ecn). In order to limit the computational cost 
of each time step, we also modify At so that we only need 
Sobj tridiagonal solves for each Crank-Nicholson time step 
(in our simulations, s bj = 14). 

To evolve the disc-magnetic star system, we also need to 
choose an implementation of the inner and outer boundary 
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conditions. We encounter two types of boundary conditions: 
Dirichlet conditions of the type y = j/bc are enforced by 



replacing the discretized version of (All by yo,N-i = Vbc, 
while Neumann conditions of the type y = y' BC are enforced 
by explicitly replacing y' by y' BC whenever necessary in ( Al I. 
If a second derivative is required to evaluate ( Al I at the 
boundary, y_i and yjv are obtained using (yi+i — yi-i) = 
(Ax)y' i and the known value of y' at the boundary. 
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